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Abstract 


Sampling  theory  is  that  branch  of  mathematics  that  seeks  to  reconstruct 
functions  from  its  values  at  a  discrete  set  of  points.  The  fundamental  result 
in  sampling  theory  known  as  "Shannon's  sampling  theorem"  has  many  ap¬ 
plications  to  signal  processing  and  communications  engineering.  I  demon¬ 
strate  Shannon's  result  via  complex  interpolation  methods.  I  then  quote  a 
result  that  uses  these  methods  to  solve  interpolation  problems  on  unions  of 
noncommensurate  lattices,  which  are  created  via  a  specific  number  of  the¬ 
oretic  guidelines.  These  interpolations  give  Shannon-type  reconstructions 
on  these  lattices.  I  close  by  doing  simulations  in  MATLAB  of  the  sampling 
reconstructions  on  these  noncommensurate  grids. 
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1.  Introduction 


Sampling  theory  is  the  distinctive  branch  of  mathematics  that  attempts 
to  solve  the  interpolation  problem  of  a  function  with  fixed  growth  from 
known  sampled  values  by  reconstructing  the  original  function.  The  the¬ 
ory  is  a  subset  of  the  general  theory  of  interpolation,  which  constructs 
functions  that  satisfy  the  known  values.  Sampling  theory  uses  interpola¬ 
tion  with  other  knowledge  of  a  function  to  find  methods  that  reconstruct 
the  original  function.  This  involves  restrictions  on  the  samples  of  known 
values. 

The  most  significant  contribution  to  this  theory  is  called  the  "classical  sam¬ 
pling  theorem."  This  theorem  is  known  by  a  myriad  of  names — Whittaker, 
Kotel'nikov,  and  Shannon— but  the  result  is  essentially  the  same.  The  theo¬ 
rem  tells  us  that  if  a  signal  /  is  a  function  of  finite  energy  on  R  (/  G  L 2  (R)) 
with  Fourier  transform  /(ta)  =  0  for  all  |w|  >  fl,  and  if  T  <  l/2f2,  then  for 
all  t  G  M, 


OO 

f(t)=T  J2  f(nT ) 

n=—o o 


sin(|?(t  —  nT)) 
it{t  —  nT) 


(1) 


In  other  words,  if  /  is  sampled  uniformly  at  a  sufficiently  high  rate,  /  can  be 
reconstructed  from  these  samples  by  the  formula  in  equation  (1).  This  for¬ 
mula  has  a  rich  history  extending  back  as  far  as  the  1600s.  However,  much 
of  the  crucial  research  in  developing  and  establishing  the  theory  has  oc¬ 
curred  in  the  twentieth  century,  primarily  because  of  the  research  in  the 
various  engineering  fields.  With  the  motivation  of  obtaining  useful  and 
practical  results,  the  classical  sampling  theorem  has  been  generalized  and 
extended  to  a  number  of  applications. 
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2.  Preliminaries  in  Harmonic  Analysis 


Before  I  discuss  the  main  results  of  this  report,  I  first  introduce  a  number 
of  preliminary  definitions  and  theorems.  Much  of  what  follows  is  standard 
notation  in  the  literature  and  can  be  found  in  any  text  covering  Fourier 
analysis.  (For  example,  see  Benedetto  [1],  Dym  and  McKean  [2],  or  Korner 
[3].  Also  see  Apostol  [4],  Higgins  [5],  Marks  [6],  and  Zayed  [7].)  First,  I 
define  in  the  following  what  are  absolutely  and  square  integrable  functions: 

Function  f.  A  function  /  is  called  absolutely  integrable,  i.e.,  /  €  X^R),  if 

[  \f(x)\  dx  =  ll/Hi  <  +oo  . 

Jr 

If  /  is  in  X1,  one  can  say  that  ||/|  |i  is  the  L 1  norm  of  /.  Similarly,  a  function 
is  called  square  integrable,  i.e.,  /  €  X2(R),  if 

[  \f{x)\2dx  =  ||/||2  <  +00  . 

J  R 

If  /  is  in  X2,  one  can  say  that  1 1/| |2  is  the  X2  norm  of  /. 

In  this  section,  all  functions  are  considered  absolutely  and  square  integrable 
functions  on  the  real  line,  unless  otherwise  noted.  Likewise,  all  integrals 
are  assumed  to  be  over  the  whole  domain  (either  R  or  C,  depending  on  the 
context)  unless  noted  otherwise.  The  Fourier  series  and  Fourier  transform 
play  integral  parts  in  the  sampling  theory  and  Fourier  analysis.  Their  def¬ 
initions,  from  Benedetto  [1]  and  Dym  and  McKean  [2],  are  defined  in  the 
following: 

Fourier  Series.  Let  /  be  a  periodic,  integrable  function  on  M,  with  period  2f l. 
Then  the  Fourier  series  of  /  is  represented  as 

OO 

/(*)=  Cnex p-wn  5  (2) 

n=—oc< 

where  the  Fourier  coefficients  cn  are  defined  by 

cn  =  ±iJQJ(x)expi™x^dx  .  (3) 

Fourier  Transform  and  Inversion  Formula.  The  Fourier  transform  of  /  €  X*(R) 
is  defined  as 


/M=  /  /(*)exp-2^di  , 
J  R 


(4) 
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and  its  inversion  formula  for  /  €  LX(R)  is 

f(t)=  [  /Mexp2™1^  .  (5) 

Jr 

The  factor  2i r  in  the  exponential  simplifies  Plancheral's  formula,  i.e.,  1 1/|  I2  = 

|  [/I |2.  Often,  the  Fourier  transform  of  /  is  also  denoted  F.  The  Fourier  trans¬ 
form  can  be  extended  to  square  integrable  functions  via  a  density  argument 
on  C^(R),  the  space  of  infinitely  differentiable,  rapidly  decreasing  func¬ 
tions  on  R.  The  transform  is  an  isometry  in  L2(R).  By  itself,  the  Fourier 
transform  has  useful  applications  in  signal  and  image  processing.  If  I  con¬ 
sider  the  function  /  to  be  a  signal,  then  I  refer  to  the  domain  of  /  as  the 
time  domain  and  the  domain  of  the  transform  as  representing  the  signal 
frequency.  The  following  definitions  are  useful  when  referring  to  the  time 
or  frequency  domain  of  a  signal  or  function: 

Support.  The  support  of  a  function  /  :  R  -»  R,  denoted  supp (/),  is  the 
closure  of  the  set  on  which  /  is  nonzero,  i.e.,  supp(/)  =  {x  £  R  :  f(x)  ^  0}. 
The  function  /  is  said  to  have  compact  support  if  supp(/)  is  a  compact  set. 

Band-limited.  The  function  /  :  R  ->  R  is  called  ft  band-limited  if  the  support 
of  its  Fourier  transform  /  is  contained  in  the  bounded  interval  [— Cl,  ft]  in  R, 
i.e.,  supp(/)  C  [-ft,  ft]- 

The  Fourier  transform  has  a  multitude  of  mathematically  convenient  and 
useful  algebraic  properties:  linearity,  symmetry,  conjugation,  translation 
(time  shifting),  modulation  (frequency  shifting),  and  time  dilation  [1].  Also, 
the  transform  has  a  number  of  important  analytic  properties:  boundedness, 
continuity,  Riemann-Lebesgue  Lemma,  time  differentiation,  and  frequency 
differentiation  [1-3,5].  One  of  the  more  relevant  properties  of  the  Fourier 
transform  involves  a  mathematical  operation  called  a  convolution. 

Convolution.  For  /,  g  €  LJ(R),  the  convolution  of  /  and  g  is  defined  by 

h(x )  =  /  *  g{x)  =  J  f(y)g{x  -  y)  dy  (6) 


for  x  G  R. 

Convolution  is  a  mathematical  operation  that  consists  of  the  filtering  of 
one  signal  (the  source)  passing  it  through  a  filter,  which  is  simply  another 
signal. 

Convolution  is  commutative,  i.e.,  f*g  =  g*f,  and  if  //  <?/  and  h  are  of  com¬ 
pact  support,  convolution  is  associative.  The  subject  of  convolution  is  im¬ 
portant  to  signal  reconstruction.  The  operation  is  the  filtering  of  one  func¬ 
tion  by  passing  it  through  a  filter,  which  is  simply  another  function.  More 
on  this  topic  will  be  discussed  in  section  4.1.  For  now,  I  note  the  following 
two  significant  properties  from  the  Fourier  perspective: 

7*9  =  /•?>  and 

7-9  =  7*9- 
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These  properties  prove  to  be  extremely  useful  in  developing  the  theory  on 
deconvolution. 

Finally,  I  conclude  with  the  Paley- Wiener  theorem: 

Paley-Wiener.  The  Fourier-Laplace  transform  of  an  infinitely  differentiable 
function  /  with  compact  support  contained  in  {|t|  <  A}  is  an  entire  func¬ 
tion  f(()  in  C,  which  satisfies  the  following  property: 

For  every  integer  N  >  0,  a  positive  constant  C  =  C(N)  exists,  such  that 

1/(01  <  £(1  +  |C|)_iVexp2^lImCi  for  all  C  €  C  .  (7) 

Conversely,  every  entire  function  in  C  satisfying  this  property  is  the 
Fourier-Laplace  transform  of  a  C°°  function,  with  compact  support  within 
{\t\  <  A}. 

The  classical  Paley-Wiener  theorem  says  that  a  square-integrable  complex¬ 
valued  function,  defined  over  the  real  line,  can  be  extended  off  the  real  line 
as  an  entire  function  of  exponential  type  if  and  only  if  its  Fourier  transform 
f(x )  is  identically  zero  for  |x  >  a,  i.e.,  if  and  only  if  /  is  band-limited  to 
[—a,  a].  (For  a  derivation,  see  literature  by  Dym  and  McKean  [2]. ) 

The  Shannon  series  is  a  means  of  extending  to  the  entire  complex  plane  C. 
The  extension  of  the  Paley-Wiener  theorem  to  generalized  functions  (to 
tempered  distributions)  is  called  the  Paley-Wiener-Schwartz  theorem.  Ad¬ 
ditional  background  material  can  be  found  in  literature  by  Benedetto  [1], 
Dym  and  McKean  [2],  and  Korner  [3]. 


3.  Classical  Sampling  Theory 


3.1  Classical  Sampling  Theorem 

Sampling  theory  is  an  area  of  high  interest  and  interacts  with  numerous 
other  mathematical  fields.  The  idea  behind  sampling  is  to  convert  continu¬ 
ous  time  signals  into  appropriate  discrete  signals  that  represent  the  original 
signal  and  vice  versa.  Reconstruction  reverses  the  process  to  retrieve  the 
original  signal  from  the  samples.  This  theory  of  sampling  was  created  via 
an  amalgam  of  mathematical,  scientific,  and  engineering  research.  Without 
dispute,  the  fundamental  result  from  the  conglomeration  of  these  sources 
has  been  what  is  called  the  classical  sampling  theorem:  Let  /  be  a  function 
of  finite  energy  on  1  (/  6  L2(M))  with  Fourier  transform  f(uj)  =  0  for  all 
|u;|  >  fi,  i.e.,  fit)  is  fi  band-limited. 

If  T  <  1  /2f I,  then  for  all  t  €  M, 


OO 

m  =  T  /(nT) 

7?  — —  CO 


sin(^(t  —  nT )) 
7r  ft  —  nT) 


(8) 


If  T  <  1  /2Q  and  /(nT)  =  0  for  all  n  e  Z,  then  /  =  0.  The  rate  ^  is  known 
as  the  Nyquist  rate. 

This  theorem,  though  it  is  simple  to  state,  has  far-reaching  consequences. 
The  theorem  basically  states  that  some  signals,  or  functions,  can  be  per¬ 
fectly  represented  by  sample  values  of  the  function  taken  at  regular,  or  eq- 
uispaced,  intervals.  Knowledge  of  the  frequency  bound  is  required,  since 
this  determines  the  minimum  rate  at  which  the  signal  must  be  sampled 
to  reconstruct  the  signal  completely.  Because  Nyquist  discovered  this  min¬ 
imum  rate,  it  is  commonly  referred  to  as  the  Nyquist  rate  [8].  No  infor¬ 
mation  is  lost  if  the  signal  is  sampled  at  this  rate,  and  also  no  additional 
information  is  acquired  by  sampling  faster  than  this  rate.  That  is,  the  signal 
is  already  perfectly  reconstructed.  Essentially,  the  band-limited  condition 
guarantees  that  sufficiently  close  sampling  reveals  a  nicely  behaved  func¬ 
tion  between  the  sample  points.  Since  /  is  band-limited,  one  can  gather 
data  sufficiently  faster  than  the  signal's  oscillation  in  time  and  then  recover 
/  using  the  "sine"  function  as  an  interpolator.  Thus,  one  can  think 

of  classical  sampling  as  a  uniqueness  theorem  from  the  theory  of  analytic 
splines. 

There  is  also  a  duality  in  the  theorem.  Not  only  can  one  reconstruct  the 
function  in  time,  assuming  the  function  is  band-limited,  but  also  in  fre¬ 
quency,  assuming  that  the  function  is  time-limited.  In  this  case,  the  function 
and  its  transform  switch  roles. 

There  is  an  intriguing  history  involving  the  development  of  this  theorem. 
A  dispute  prevailed  over  who  initially  discovered  this  result  and  who  is 
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credited.  Consequently,  the  theorem  has  several  names  associated  with  it: 
the  cardinal  series,  the  Whittaker  sampling  theorem,  the  Kotel'nikov  theo¬ 
rem,  and  the  Shannon  sampling  theorem.  Often,  the  names  are  combined 
in  various  forms.  For  simplicity,  it  is  referred  to  in  this  report  simply  as  the 
classical  sampling  theorem. 


3.2  A  Proof 


There  are  numerous  proofs  of  the  classical  sampling  theorem  [1,5-7,9,10], 
some  more  rigorous  and  satisfying  than  others.  A  simple  proof  of  the  theo¬ 
rem  follows: 

Proof.  Let  /  be  O  band-limited,  i.e.,  supp(/)  c  [— 0,  O].  Since  this  condition 
holds,  /  can  also  be  represented  as  the  restriction  of  a  20  periodic  function 
to  the  support  of  /  as 

OO 

/(«)=  Cn exV~mnu/il  'x{-n.n\iu)  >  (9) 

71—  —  OO 


where  the  Fourier  coefficients  are  simplified  with  the  use  of  the  Fourier 
inversion  formula  and  the  nature  of  X.  Note,  Xa  is  the  characteristic  func¬ 
tion  defined  to  be  unity  in  the  set  A  and  zero  elsewhere.  Then  the  Fourier 
coefficients  are  found  to  be 


Cn  =  T  y" /(u)ex du 

=  - -j  /(«)  exp2^(n/2n)u  du 

(10) 

— /(— )  . 
20jy20’ 


Substituting  equation  (10)  back  into  equation  (9)  and  solving  for  /  with 
the  use  of  the  inverse  Fourier  transform  in  equation  (5)  and  the  uniform 
convergence  of  Fourier  series,  the  classical  sampling  theorem  is  the  result 


/(*) 


(uj  exp 


/OO 

/(t 

-OO 

/OO  1  00  ,n 

255  E  /(^)exP"”n“/n^[-aniexp2’i“’<i« 

■°°  n=— oo 

OO  rQ 

25)  E  A| 5>  /  ex p<2— 


(11) 


n=— oo 
00 


-1  y  /(il)2!!-5in,rl2!!l-’l> 

2 O  ^ 


E  /W 


2  JV  7r(2  Q,x  —  n) 

n  sin7r(2flx  —  n) 

2fl  7t(2Qx  —  n ) 


The  third  equality,  where  the  integral  and  summation  are  exchanged,  is 
permitted  by  uniform  convergence.  So,  if  T  =  1/2Q,  then  equation  (11)  is 
identical  to  equation  (1). 
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3.3  An  Example 


The  sampling  theory  with  this  band-limited  condition,  at  first,  appears  too 
restrictive  for  practical  use.  The  formula 


OG 

/(f)  =  T  ^2  f(nT) 

n= — oo 


sin(^(f  —  nT )) 
ir(t  -  nT) 


requires  that  the  signal's  frequency  be  bounded  for  the  reconstruction  to  be 
stable.  However,  in  practice,  this  condition  is  not  restrictive.  The  physical 
interpretation  of  an  Q  band-limited  function  is  that  it  is  a  signal  with  a 
frequency  no  greater  than  f!  cycles  per  second.  In  countless  examples,  this 
reflects  the  reality  of  the  physical  world.  For  instance,  many  sounds  are 
band-limited  naturally;  i.e.,  the  frequency  of  orchestral  music  is  less  than 
20  kH  and  the  human  voice  is  less  than  8  kH.  Even  our  hearing  is  limited 
to  a  certain  range  of  frequencies,  and  outside  of  these  frequencies,  we  miss 
the  signal.  For  instance,  not  many  people  can  hear  the  frequency  of  a  dog 
whistle. 


Suppose,  one  had  a  signal  that  was  f!  band-limited.  Then  its  Fourier  trans¬ 
form  would  be  equivalent  to  the  product  of  S(u>)  and  F(u)  =  ^(w), 

the  characteristic  function  of  [— Then  the  signal  is  defined  by  the 
convolution  of  the  inverse  Fourier  transform  of  S(oj),  perhaps  s(x),  and 
f(x)  =  ^|— .  For  simplicity,  let  Cl  =  1  and  the  function  S(uj)  =  1,  then  the 
signal  is  represented  just  by  f(x).  (See  fig.  1.) 

According  to  the  sampling  theorem,  since  /  is  band-limited,  it  can  be  re¬ 
constructed  from  its  values  at  the  points  nT,  where  T  <l/2fI  =  l/2  and 

71  G 


oo 

ff(x)  =  T 

n=— oo 


f(nT ) 


sin(^(x  —  nT)) 
tt(x  —  nT) 


where  x  €  R. 

In  the  reconstruction  in  figure  2, 1  set  T  =  l/2fi  =  1/2.  Because  of  the  band- 
limited  condition,  there  is  no  gain  to  oversampling.  To  further  illustrate 
the  significance  of  the  band-limited  condition,  i.e.,  the  need  to  sample  at 
or  above  the  Nyquist  rate,  the  signal  function  is  under  sampled.  The  next 
sample,  in  figure  3,  is  taken  at  T  =  be.,  below  Nyquist. 

Knowledge  of  the  frequency  bound  allows  one  to  periodically  extend  the 
transform  to  the  real  line.  But  when  the  sample  is  taken  below  the  Nyquist 


Figure  1.  Original  signal 
f(x)  and  its  Fourier 
transform 
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Figure  2.  Original  signal 
/  and  reconstruction  g. 
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Figure  3.  Aliasing 
caused  by 
undersampling. 
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rate,  those  periodic  extensions  contain  overlapping  sections  and  the  infor¬ 
mation  is  not  weighted  appropriately.  This  effect,  demonstrated  in  figure  3, 
is  known  as  aliasing. 

3.4  Deriving  the  Cardinal  Series 

There  are  numerous  ways  in  which  to  derive  the  cardinal  series,  or  sam¬ 
pling  formula,  as  attested  in  the  literature  concerning  this.  The  histori¬ 
cal  development  of  its  construction  follows  in  section  3.5.  In  this  report, 
I  display  several  methods  that  successfully  derive  the  cardinal  series.  The 
first,  which  is  called  the  delta  method  [5],  is  natural  in  its  approach  but  is 
not  rigorous!  The  second  derivation  uses  the  Cauchy  integral  formula  and 
shows  this  to  be  equivalent  to  the  classical  sampling  theorem  under  certain 
conditions. 


3.4.1  Delta  Method 

For  this  derivation,  the  function  is  sampled  equidistantly.  Mathematically, 
this  is  using  a  series  of  generalized  functions  with  point  support,  i.e.,  delta 
functions.  Thus,  a  sampled  version  of  /  is 

fs{t)  =  /(t)^5(t-nr)  =  ^/(nr)6(f-nr)  ,  (12) 

n&  n€Z 

where  r  controls  the  rate  of  the  sampling.  This  hints  at  the  following  recon¬ 
struction  formula  that  uses  some  function  g: 


m  =  E  f(nT)g(t-nT )  .  (13) 

neZ 

The  task  is  now  to  find  an  appropriate  interpolating  function  g  and  rate  r . 
Then  equation  (13)  can  be  rewritten  as 


f(t)  =  f  f{u)g(t  -u)^2  S(u  -  nr)  du  .  (14) 

neZ 

The  Fourier  series  expansion  of  the  delta  sum  is  Ylne ~  nr)  = 
1/r  Yju&l  exp_2iri7lu/T  with  the  use  of  equation  (12)  [1].  So,  again  equation 
(14)  is  rewritten  as 


m 


[  f(u)g(t-u)1-'£/exp-2'linu/Tdu 
^  T  n£7. 

/(•)X>  P-2,rin/T*^(-)(*). 

n€Z 


(15) 


where  the  convolution  operator  *  is  defined  in  section  2.  Then,  the  Fourier 
transform  of  both  sides  results  in 

J{x)  =  i<?(x)  f{x  +  •  (16) 

T  neZ 
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Thus,  the  transform  of  /  is  equivalent  to  some  weighted  spectral  repeti¬ 
tion  of  itself.  To  make  sense,  g  must  be  such  that  it  chooses  only  one  copy 
of  the  spectrum.  This  is  inherently  impossible  if  the  repetitions  overlap. 
Supposing  that  /  is  u;  band-limited,  then  \/r  >  2cc  is  needed  to  prevent 
overlapping  repetitions.  This  reveals  the  importance  of  the  rate  r  <  1/2 to 
in  sampling  the  function.  Now,  g(x)/r  is  a  window  in  which  to  keep  the 
desired  copy  of  the  spectrum,  where  n  =  0.  Thus,  define  g  so  that 


. L](«)  ,  (17) 

and  note  that  g  needs  only  to  be  equivalent  to  the  characteristic  function 
over  the  interval  [— oo,u]  c  [— l/2r,  1/2t].  This  is  equivalent  to  equation 
(17)  and  keeps  the  consistency  in  the  formula  dependent  on  the  sampling 
rate  exclusively.  This  reveals  the  interpolator  to  be 


g(t)  =  t 


sin  2-K-^t 
irt 


sin  jt 

H  ’ 


(18) 


and  therefore,  the  final  rewriting  of  equation  (13)  becomes  the  sampling 
formula 


f{t)  = 

n€Z 


=  X^(nr) 


sin  f(t  —  nr) 
f(t-  nr) 


(19) 


Again,  the  classical  sampling  formula  is  derived.  This  is  not  the  only 
method  of  deriving  the  classical  sampling  theorem.  A  method  that  uses 
complex  variables  is  demonstrated  in  the  next  section  (3.4.2). 


3.4.2  Complex  Method 

The  sampling  theorem  is  intimately  connected  to  a  wide  variety  of  mathe¬ 
matical  theories.  Under  certain  conditions,  it  has  been  shown  to  be  equiva¬ 
lent  to  a  number  of  other  important  formulas:  the  Poisson  summation  for¬ 
mula,  the  Cauchy  integral  formula,  etc  [5].  Its  connection  to  Cauchy  inte¬ 
gral  theory  is  demonstrated  here. 

The  relation  between  complex  integration  theory  and  the  sampling  formula 
is  inherently  obvious.  The  contour  integral  around  a  set  of  simple  poles 
is  simply  a  linear  combination  of  some  function's  values  at  the  poles  (of 
course,  with  the  poles  removed).  This  is  a  kind  of  sampling  relationship, 
made  obvious  when  one  considers  that  in  the  classical  sampling  theorem, 
the  "sine"  function  contains  the  appropriate  zero  values  at  the  sampled 
data  points  of  the  signal  function.  This  technique  fits  under  the  umbrella 
of  general  complex  interpolation.  The  splitting  of  the  integral  is  called  the 
Jacobi  interpolation  formula.  This  technique  was  used  by  Berenstein,  Casey, 
Yger,  and  Walnut  to  solve  deconvolution  problems. 

To  show  this  connection  mathematically,  I  will  start  with  Cauchy's  integral 
formula.  See  Marsden  and  Hoffman  [11]  for  the  necessary  background  in 
complex  analysis. 
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Cauchy  Integral  Formula  Theorem  [11].  Let  /  be  analytic  on  a  region  D  and  let 
C  be  a  closed  curve  in  D  that  is  homotopic  to  a  point,  i.e.,  let  C  be  a  simple, 
closed,  and  rectifiable  curve  in  C  and  also  let  z  €  D  be  a  point  in  D  inside 
C.  Then,  for  r  €  No, 

,)(c,  2)  ■  /wW  =  £-fc  (Triym  dl  ■  <20> 

where  rj(C,  z)  =  l/2m  §c  \/{y-z)  dy  is  the  winding  number.  Note,  rj{C,  z)  = 
1  if  C  is  positively  oriented  and  loops  around  2  €  int  C  only  once. 

Given  this  theorem,  I  can  derive  the  classical  sampling  formula.  For  sim¬ 
plicity,  assume  that  Q  <  1/2 T  and  place  the  sample  values  at  the  integers, 
n  g  Z.  Thus,  the  sampling  rate  T  =  1  and  <  1/2.  Let  /  be  a  function 
that  satisfies  the  conditions  in  the  theorem.  By  the  Paley-Weiner  theorem 
in  section  2,  /  is  an  analytic  function  int  6  R,  which  analytically  contin¬ 
ues  to  the  entire  complex  plane.  Moreover,  this  entire  function  satisfies  the 
Paley-Weiner  growth  bound,  i.e.,  for  Fl  >  0, 

.  (21) 

Now,  define  /2(C)  =  sin/rC.  This  has  zeros  at  the  integer  values,  i.e.,  the 
sampling  rate.  So  the  contour  will  need  to  avoid  the  zeros  of  //(C)-  I  define 
Cm  to  be  a.  square  contour  with  comers  on 

(±l±t)(m+  i)  • 

Thus,  r)(C,  z )  =  1.  Now,  first  use  the  Jacobi  interpolation  formula  (see  Casey 
and  Walnut  [12])  that  gives 


m  = 


*  JUr, 


(1)  + 


l 

2iti 

1 

Jcm  C 


(2) 


Next,  I  shall  show  that  (2)  -^Oasm 

y{z) 


(2)  = 


Jc, 


00,  where 

/(C) 


d( . 


(23) 


2vr*  JCm  (C  -  *)£«) 

Recall,  that  Cm  is  a  square  contour  with  corners  (±1  ±i)(m  + 1/2).  Let  R  be 
a  closed,  compact  disc  centered  at  0.  Then  when  m  +  1/2  >  R  and  C  €  Cm, 
|C  —  z\  >  m  +  1/2  -  R,  such  that 

,m(z)  I  /(C) 


1(2)1 


Ja 


< 


< 


27li  Jcm  (C  -  z)KO 
IPWI  /  1/(01 


d0 


27T 


£ 

l/SWI 


\(c-z)\-\m\ 


27r(m  + 


Li 

b  -  R)  Jc, 


1/(01 

Im(C)I 


d( 


dC- 


(24) 
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Now,  let  C  =  u  +  iv.  Then,  define  Cm  as  the  side  of  Cm  that  joins  (1  — 
i)(m  +  1/2)  and  (1  +  i)(m  +  1/2),  Cm  that  joins  (1  +  i)(m  +  1/2)  and  (-1  + 
i)(m  +  1/2),  Cm  that  joins  (-1  +  i)(m  +  1/2)  and  (-1  -  i)(m  +  1/2),  and 

Cm  that  completes  the  square  contour  that  joins  (— 1  —  i)(m  +  1/2)  and 
(1  -  i)(m  +  1/2). 

On  the  side  Cm\  (  =  m  +  1/2  +  iv ,  thus 

|  sin(7r<^)|2  =  sin2(7ru)  +  sinh2(7rv)  =  1  +  sinh2(7ri>)  =  cosh2(7rv)  >  exp2?r^  . 


Also  |/(C)|  <  A  exp  27rQH  in  equation  (21),  so  that  if  Q  < 

<f  ¥MdC  <  A  [m+  2  dv 

Jett  KOI  J-m-\ 

exp_7r(2n+1^m+^  -1 
=  A  tt(2Q  +  1)  ’ 

which  disappears  as  m  — >  oc,  since  Q  >  0.  However,  if  Q  =  1/2,  then 

m+i 


(25) 


i  i/«)i ,,  „  r+5  m 
Jcs>  im(oi  c  -  /-„-i  a+icir  c 


=i: 


(1  +  (m  +  i)2  +  v2)N 


dv 


(26) 


A\  ■  (2m  +  1) 


(1  +  (m  +  \)2)N 

Thus,  the  integral  goes  to  0  as  m  — +  oo.  A  similar  approach  can  show 
the  same  for  C(m}.  On  the  side  C$ ,  C  =  u  +  i(m  +  1/2),  and  |  sin(7rC)|  > 
sinh(7r(m  +  1/2)).  Thus,  if  ft  <  1/2,  then 

^-27r\m+~\Ct 


i  <  r+* 

Jett  KOI  J-m- i 


A! 


exp 


(1  +  u2  +  (m  +  ^)2)N  sinh  tt (m  +  |) 

A\  •  (2m  +  1)  exp~27rn(m+f) 

(1  +  (m  +  l)2)^  sinh7r(m  +  5) 


du 


(27) 


which  also  disappears  as  m  —>  00,  since  fi  >  0.  Likewise,  Cm1  is  integrated 
in  the  same  manner.  So,  equation  (22)  becomes 


/(*)  =  (1)  = 


/( 0  £(0  -  Kz) 


Jcm(~z 

Kz) 

27 n 


M(C) 


dC 


/  /«) . 

l 

1 

tcm  C-z 

ji(0 

M(0. 

dC. 


At  this  point,  I  need  to  note  a  couple  items  before  proceeding.  First, 


N 


sin7rC  =  lim  s/v(0  =  lim  nC  TT  f  1  +  - 

N—+OO  /V— ™  1  n 


N—+00 


(28) 


(29) 


3= 1 
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is  the  product  representation  of  a  sine  function.  Second,  note  the  following 
partial  fractions  decomposition: 


1 

(C-  z)sN{() 


1 

(C  -  z)sN{z) 


+  |nH  (C  ~  n)(n  "  2)S'w(n) 


(30) 


To  show  this,  I  need  the  following:  Let  p(()  be  a  polynomial  with  only  sim¬ 
ple  zeros  at  {ry}j=1,  then  (shown  in  app  A) 


1 

p(0 


E 


3= 1 


1 

(C  -  rj)p'(rj) 


(31) 


Then,  if  m  <  N  and  defining  Rrn  =  (2)  on  Cm, 


m  = 


2tt  i 


/(C)  . 

1 

1 

tm  C  -  * 

jl{z) 

m(C)_ 

=  i  m 


-f 

'*  JCr 


N—>oo  2m  Jcm  C 

£  E 


dC  +  Rrn 

1 _ 1_ 

Sjv(^)  Sjv(C) 

/(C) 

iV-io  2tt*  (C-n)(n-z)s)v(n) 


=  Um  SN<2) 


=  lim 


sN{z 


p  E  ■/ 

1  ,-T^r  JCr 


/(C) 


d(  +  Rm 

dQ  +  Rrn 

dC  +  Rm  (32) 


—  lim  sn(z)  / 

AT-^oo  ' 


/(n) 


—  Sill  TCZ 


\n\<m 


(z  -  n)s’N(n ) 


+  Rn 


( _ I  \n 

V  /(n)- 
^  7r(>6  —  n) 


|n|<m 

=  ^  /(n)sinc(z  —  n)  +  . 

|n|<m 


Thus,  as  m  -»  oo,  then  -*  0,  and  this  is  the  classical  sampling  theo¬ 
rem.  After  the  limit,  the  contour  integral  covers  over  all  the  zeros  of  the 
sampling  function  and  thus  integrates  over  the  whole  function.  Note,  the 
switch  between  limit  and  integration  in  the  second  equality  and  between 
summation  and  integration  in  the  fourth  is  permitted  because  of  uniform 
convergence.  A  similar  approach  is  undertaken  by  Higgins  [5],  and  he  also 
demonstrates  other  and  more  general  equivalence  relations  involving  the 
classical  sampling  formula. 


3.5  A  Bit  of  History 

Sampling  theory  grew  out  of  a  plethora  of  other  mathematical  fields.  It  is 
intimately  connected  with  interpolation  theory  and  approximation  theory. 
In  fact,  for  certain  cases,  it  can  be  shown  that  the  cardinal  series  (the  classi¬ 
cal  sampling  formula),  the  Cauchy  integral  formula,  and  the  Poisson  sum¬ 
mation  formula  are  equivalent  [5].  Sampling  theory  also  has  applications 
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in  numerous  other  fields  of  systems  theory,  prediction  theory,  information 
theory,  stochastic  processes,  signal  and  image  processing,  etc.  The  theory 
of  interpolation  grew  out  of  the  need  to  calculate  intermediate  values  of 
functions  with  some  known  values.  In  1624,  Briggs  introduced  methods 
for  achieving  this  by  using  successive  and  modified  differences.  Wallis  also 
had  some  influence  in  this  area  in  the  1650s  and  is  credited  with  coining 
the  term  "interpolation."  However,  it  was  not  until  the  work  of  Gregory 
and  Newton  that  the  polynomial  nature  of  these  methods  began  to  be  rec¬ 
ognized.  In  1670,  Gregory  introduced  the  first  interpolation  series,  and  later 
Newton  introduced  divided  and  adjusted  divided  differences.  Both  meth¬ 
ods  derive  from  taking  the  polynomial  as  a  close  approximation  to  the  (con¬ 
tinuous)  function  to  be  interpolated  [13]. 

Although  these  methods  were  somewhat  successful  in  providing  approxi¬ 
mations,  it  became  desirable  to  construct  a  polynomial  interpolant  without 
the  need  of  finding  successive  differences.  This  was  achieved  with  the  iden¬ 
tity  below: 


/(Ao)  +  (t  —  Ao)/(Aq,  Ai)  + - 1-  (t  —  Aq)  . . .  (t  —  Am-i)/(Ao,  . . .  ,  Am) 


M 


=£  m 

i= o 


Gm( t) 


(33) 


where 


/( Ao.Ai) 

/(Aq,  Ai,  A2) 


/(Ai)-/(Aq) 

Ai  —  Ao 

/(Ai,A2)  -  /(Aq,  Ai) 
A2  —  Aq 


and  . 

M  . 

aM(t)  =  ft  (1  -  j- 

j— 0  v  3 

Both  sides  are  identical  because  both  are  polynomials  of  degree  M  taking 
on  the  same  values  on  (A j}.  The  right  side  of  the  equality  in  equation  (33) 
is  known  as  Lagrange's  formula.  Essentially,  the  Lagrange  interpolator  acts 
as  a  delta  function  over  the  known  values.  This  formula  first  appeared  in 
Lagrange's  lectures  given  in  1795.  However,  it  was  previously  discovered 
by  Waring  in  1779  [5]. 

Parallel  to  the  above  approach,  it  became  desirable  to  have  a  periodic  inter¬ 
polant  rather  that  a  polynomial  one.  Given  a  periodic  function  of  the  form 
p{t)  =  Z]L_„  Cj  expyt,  then  appropriate  formulas  were  found  for  these 
cases.  One  such  formula  is 


p(t)  = 


2A/+1 


_ )(t _ 

2  )\L  2M+1 


)] 


gill  —( /  —  _ ^ 

bm  2  2M+1 1 


2M+  1 


This  finite  sampling  series  is  credited  to  Cauchy  (1841),  although  Gauss  in¬ 
troduced  primitive  forms  as  far  back  as  1805  [5].  Cauchy's  other  significant 
contribution  was  to  be  the  first  to  note  the  importance  of  the  rate  of  the 
sampling. 

The  research  conducted  in  the  nineteenth  century  extended  Lagrange's  for¬ 
mula  from  interpolation  on  a  finite  number  of  points  to  the  more  general 
case  of  infinitely  many  points  {An}.  Assuming  / (An)  =  an  for  all  n  e  N, 

ff\_  </l(-g) 

.  /^-^GV(A„)(,-An)  1 

where  <f>  is  a  function  with  simple  zeros  at  {A„}.  A  more  general  formula 
was  also  constructed: 

n  (  —  )sn _ - 

An  <t>'(\n)(z-\n)  ’ 

where  (s„)  is  a  sequence  of  integers  appropriately  chosen  to  ensure  conver¬ 
gence  [14].  These  formulas  were  found  useful  by  Cazzaniga  in  1882  for  con¬ 
structive  purposes.  The  century  ended  with  the  first  explicit  statement  of 
the  cardinal  series  by  Borel.  In  1899,  Borel  set  An  =  n  and  noted  that  sin(Trz) 
has  the  appropriate  zeros,  returning  the  classical  sampling  formula  [5]. 

During  the  first  half  of  the  twentieth  century,  this  special  case  of  the  La¬ 
grange  formula  with  knowledge  of  equidistantly  spaced  points,  i.e.,  the 
cardinal  series,  was  rediscovered  in  this  form  several  times  [5,14].  Whipple 
did  so  in  1910,  and  E.  T.  Whittaker  in  1915  [10].  Whittaker  would  also  point 
out  the  band-limited  nature  of  the  sum.  His  son,  J.  M.  Whittaker,  would 
coin  the  phrase  "cardinal  series"  in  1929.  Ogura,  in  1920,  was  the  first  to 
provide  a  fully  rigorous  proof  of  the  cardinal  series  by  using  the  calculus  of 
residues.  In  1928,  Nyquist  established  the  time-bandwidth  component  of 
a  signal,  the  importance  of  rates,  in  connection  with  telegraphy.  His  work 
showed  that  stability  and  reconstruction  were  not  possible  below  a  cer¬ 
tain  sampling  rate.  In  1941,  Hardy  noticed  that  the  cardinal  series  is  an 
orthogonal  expansion,  an  important  development  contributing  to  Hilbert 
and  Banach  spaces.  In  1933,  Kotel'nikov  introduced  sampling  theory  into 
communication  theory.  However,  this  fact  was  unknown  in  the  west  until 
the  1960s.  For  the  Western  world.  Shannon  introduced  it  in  1949  [14].  Shan¬ 
non's  original  statement  of  the  theorem  is  "If  a  function  / (t)  contains  no 
frequencies  higher  than  W  cps,  it  is  completely  determined  by  giving  its 
ordinates  at  a  series  of  points  spaced  (1/2  W)  s  apart." 

In  fact,  he  claims  that  this  fact  is  common  knowledge  in  communication 
and  credits  Whittaker  for  an  earlier  form  of  the  theorem.  He  also  notes 
Nyquist's  and  Gabor's  use  of  this  fact.  However,  his  name  remains  attached 
to  the  theorem. 

Since  Shannon's  use  of  the  sampling  theorem,  much  has  been  developed  in 
this  field  of  research.  Parzen  [15]  in  1956  and  Petersen  and  Middleton  [16]  in 
1962  extended  the  sampling  theorem  to  multiple  dimensions,  and  Kramer 
developed  a  generalization  of  the  classical  sampling  theorem  in  1959  [17]. 
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4.  Multirate  Sampling  Theory 


The  classical  sampling  theorem  shows  that  one  can  effectively  recover  a  sig¬ 
nal  using  a  regular  sampling  rate  determined  by  knowledge  of  the  signal 
frequency's  bandwidth.  However,  if  the  sampling  rate  is  slower  than  the 
Nyquist  rate,  the  function  cannot  be  reconstructed  via  the  classical  theo¬ 
rem.  Research  in  this  direction  has  proven  effective  via  other  methods.  The 
procedure  here  relies  on  both  the  classical  theorem  and  deconvolution  the¬ 
ory  to  show  how  a  band-limited  function  that  is  sampled  below  Nyquist 
can  still  be  completely  recovered. 

4.1  Multichannel  Deconvolution 

Deconvolution  has  proven  itself  to  be  a  useful  mathematical  tool  in  the  field 
of  signal  and  image  processing.  For  images,  it  acts  as  an  enhancing  filter  to 
correct  blurs  in  a  picture.  For  signals,  it  can  be  used  to  correct  distorted 
line  shapes  without  a  loss  of  signal-noise  ratio.  In  essence,  deconvolution 
is  useful  in  signal  and  image  processing,  where  the  incoming  data  contains 
a  high  degree  of  information.  The  reason  this  works  so  naturally  is  because 
the  convolution  equation  models  a  number  of  linear  systems. 

Before  I  progress  further,  it  is  appropriate  to  understand  what  convolu¬ 
tion  is  defined  as  mathematically.  Recall  the  definition  of  convolution,  also 
given  earlier  (sect.  2),  for  f,g  6  Lx  (R),  the  convolution  of  /  and  g  is  defined 
by 


/ 


*  sfa)  =  J  f(y)g(x  -  y) 


dy 


for  x  €  R. 

For  the  purpose  of  this  report,  the  convolution  equation  s  =  f  *  g  mod¬ 
els  linear,  translation  invariant  systems  (e.g.,  sensors,  linear  filters).  In  this 
model,  /  is  the  input  signal  function,  //  is  the  system  impulse  response  dis¬ 
tribution,  and  s  is  the  output  (received)  signal.  However,  in  many  physical 
applications,  s  is  often  a  poor  approximation  of  the  signal  /.  This  moti¬ 
vates  one  to  deconvolve  /  from  g  to  attain  the  original  signal.  Results  have 
shown  that  if  the  convolver  g  is  time-limited  (i.e.,  compactly  supported) 
and  nonsingular  (i.e.,  not  a  delta  generalized  function),  then  this  problem 
is  ill-posed,  in  the  sense  of  Hadamard  [12].  It  has  been  shown  to  be  ill-posed 
for  all  realizable  convolvers — all  convolvers  that  can  be  built.  For  the  cir¬ 
cumvention  of  this  scenario,  a  theory  of  multichannel  deconvolution  has 
been  developed  to  solve  these  equations.  A  multichannel  system  preserves 
information  about  the  signal  that  would  otherwise  be  lost.  Thus,  data  lost 
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by  one  convolver  can  still  be  retained  by  another  convolver.  The  signal  is 
now  overdetermined,  i.e.. 


Si  =  f  *  lM,i  =  1  •  •  -n  . 

Then,  if  the  convolvers  {m}  satisfy  the  condition  of  being  strongly  coprime, 
then  deconvolving  for  /  is  now  a  well-posed  problem. 

Strongly  coprime.  Strongly  coprime  is  defined  as  a  set  of  convolvers  {pi}  that 
satisfy  the  inequality 

(X>KOI2)  2  >  ^exp-BlImzl  -(1  +  Id)-  (34) 

for  every  (  e  C,  where  A  and  B  are  positive  constants  and  N  is  a  positive 
integer,  is  said  to  be  strongly  coprime. 

Note,  this  tells  us  that  l/£"=i  \MQ\2  satisfies  the  Paley-Weiner  growth 
bound.  If  this  is  so,  then  there  exists  a  set  of  time-limited  deconvolvers  {I'i} 
such  that 

Pl  *  V\  +  •  •  •  +  pn  *  vn  —  &  i 

and  consequently, 

pl  '  V\  +  •  •  •  +  pn  '  vn  =  1  i 

where  5  is  the  Dirac  delta  function.  The  second  of  these  two  equations  is  the 
analytic  Bezout  equation.  The  existence  of  such  deconvolvers  is  guaranteed 
by  the  following  theorem: 

Hormander  [12].  For  compactly  supported  distributions  {pi}f=i  on  M,  com¬ 
pactly  supported  distributions  exist  such  that 

S  =  pi  *  V\  + - 1-  pn  *  vn  . 

That  is,  they  satisfy  the  analytic  Bezout  equation  if  and  only  if  the  set  of 
distributions  {M?=1  is  strongly  coprime. 

By  Hormander 's  theorem,  a  strongly  coprime  set  has  a  solution  in  the  ana¬ 
lytic  Bezout  equation.  Thus,  given  the  deconvolvers  and  output  signals,  / 
is  naturally  produced: 


s£jsi*vi  =  *  Mi)  *  Vi  =  Y  f  *  (w  *  "0  =  /  *  *Vi)-S*5  =  f  ■ 

i  i  i  i 

The  system  is  such  that  no  information  is  lost  in  this  process.  This  occurs 
because  the  condition  of  being  strongly  coprime  guarantees  that  the  zeros 
in  the  analytic  Bezout  equation  do  not  cluster  quickly  as  |C|  ->  oo.  If  the 
{pi}  did  have  a  common  zero,  then  s'i(Q  =  0  at  that  zero  and  information 
about  /  would  be  lost.  Because  the  system  is  engineered  toward  eliminat¬ 
ing  any  common  zero,  no  information  about  /  is  lost  (and  the  problem  is 
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well-posed).  Thus,  the  signal  /  is  gathered  by  this  strongly  coprime  system, 
and  these  received  signals  are  then  filtered  by  the  deconvolvers  to  recon¬ 
struct  /. 

These  methods  are  linear  and  realizable;  thus  deconvolution  at  a  time  sam¬ 
ple  only  depends  on  the  information  near  that  time  sample.  Unfortunately, 
Hormander's  theorem  is  an  existence  theorem;  it  does  not  reveal  what  the 
deconvolvers  might  be.  One  obvious  solution  is 

_  Tm 
*  EiM2  ‘ 

The  deconvolvers  have  been  shown  not  to  be  unique,  and  so,  for  certain 
scenarios,  one  set  of  deconvolvers  may  be  better  than  another.  More  in  this 
direction  is  given  by  Casey  and  Walnut  [18]. 


4.2  Noncommensurate  Sampling  Lattices 

The  theory  developed  here  merges  together  the  two  ideas  of  multichan¬ 
nel  deconvolution  and  classical  sampling  theory.  It  has  been  developed  by 
Casey  [19],  Casey  and  Sadler  [20],  and  Casey  and  Walnut  [12,18]. 

By  the  results  in  multichannel  deconvolution,  a  set  of  strongly  coprime  con¬ 
volvers  {/-/,}  needs  to  be  created  for  the  problem  to  be  well-posed.  First,  a 
definition  is  needed,  then  the  following  theorem  will  be  useful  in  creating 
this  set. 

Definition.  A  real  number  a  is  said  to  be  poorly  approximated  by  rationals 
if  an  integer  N  >2  and  a  constant  C  —  C(a)  exist  such  that  for  all  integers 
p,  q  with  q  >  2, 


(35) 


This  class  of  numbers  will  be  denoted  by  P. 

Theorem.  Let  0  <  r\  <  •••  <  rm,  m  >  d  =  1  satisfy  the  condition  that  for 
all  i  ^  j,  ri/rj  are  in  P,  then  {X[-run]d}  is  a  strongly  coprime  set.  Then  to 
create  a  set  of  appropriate  convolvers,  one  needs  only  to  find  such  a  set 
whose  ratios  are  inadequately  approximated  by  rationals. 


4.2.1  Two  Sampling  Lattices 

Suppose  /  €  L2(R).  Let  a  be  an  irrational  in  the  class  P.  Now,  the  pair  p\  = 
p2  =  X[-a,a]  are  strongly  coprime.  Thus,  the  problem  of  solving  the 
analytic  Bezout  equation  for  v\  and  V2 

Pi  ■  v\  +  P2  ■  =  1 

is  well-posed.  Now,  pi(io)  =  sin|r2JCt'‘>  and  p2{u)  =  sin^^  have  zeros 

=  (?}  •  -  {£} 
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for  n  <S  Z\{0}.  Thus,  one  solution  to  the  analytic  Bezout  equation  is  v\  (w)  = 
l/2jJi(w)  for  w  €  Zfa  and,  likewise,  v2  (tv)  =  I/2/I2M  forw  e  Z^.  Then,  the 
problem  simply  becomes  an  interpolation  problem  on  the  set  of  zeros  T  = 
Zn ,  II  Zn„.  So,  one  has  a  new  reconstruction  formula  for  /  from  its  values 
on  T.  However,  this  can  be  done  for  any  poorly  approximated  rational.  The 
following  theorem  results: 

Theorem.  Let 

r={f}u{g}  <36> 

for  n  e  N,  and  let  A*  e  T.  Let  /  be  a  (1  +  a) -band-limited  function,  then  / 
can  be  conditionally  reconstructed  by  the  formula 

/(<)  =  £  /(x)G,(A°((f„Ai)  +  + s' mm) .  (37) 

where 


G(t )  =  sin(27rf)  •  sin(27 rat)  (38) 

km  =  Swr,  (39) 

2!  ^ 

**(*>  =  m  •  <«> 

2!  1 

Note  that  the  interpolators  at  the  origin  appropriately  model  the  function  at 
the  origin,  i.e.,  A’i  (0)  =  1  while  K[  (0)  =  0  and  AA(0)  =  0  while  K'2(0)  =  1. 
How  these  interpolators  are  calculated  will  be  demonstrated  later  in  sec¬ 
tion  4.3.1.  Also,  note  that  the  information  in  the  signal  can  be  reconstructed 
uniquely  by  its  samples  on  TU{0}.  Here,  it  is  also  important  to  note  that 
the  sampling  rates  correspond  to  1 -band-limited  functions  and  o-band- 
limited  functions.  These  add  up  to  the  band-limit  of  the  function  to  be  sam¬ 
pled.  Thus,  a  (1  +q) -band-limited  function  is  reconstructed  via  significantly 
lower  sampling  rates. 

4.2.2  An  Arbitrary  Number  of  Lattices 

This  result  generalizes  for  an  arbitrary  number  of  sampling  rates.  Let 
{rj}”=1  be  rates  such  that  all  r^/rj,  i  ^  j,  are  in  P,  the  class  of  irra¬ 
tionals  poorly  approximated  by  rationals.  Then,  {r*j}  is  a  strongly  coprime 
set.  Thus,  our  convolvers  have  the  form  m(t)  =  These  {/x,} 

model  the  impulse  response  of  a  multichannel  system.  The  convolvers  have 
Fourier  transforms 

-  ( /-\  sin(27rrjC) 

»(0  =  — 

with  zero  sets  Z*  =  {|jr},  both  respectively,  for  k  e  N.  Note  that  exclusive 
of  the  origin,  the  zeros  sets  are  nonrepetitive.  Let  =  Z,  and  T  =  (J"=i  — 
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{Aj}.  Thus,  /  can  be  reconstructed  on  its  values  at  T.  This  can  also  be  done 
via  techniques  in  complex  interpolation  theory. 

Thus,  the  following  theorem  results:  Let  {ri}"=1  be  a  set  whose  ratios  are 
poorly  approximated  by  rationals,  and  let  /  be  a  Q>2  r,) -band -limited  func¬ 
tion.  Let 


for  k  €  N  and  i  =  1, . . .  ,  n.  Also,  let 

n 

U  Ti  =  r  =  {AO  . 

z—1 

Then  /  is  uniquely  determined  by 

{/(x)}  U  {/(o),--.,/(”-!,(°)}  ■ 

Furthermore,  /  can  be  reconstructed  from  its  values  on  T  [J{0}  by  the  fol¬ 
lowing  formula: 


m  =»  E  /ft) 

A<er 


<2(0 

-  A i) 


i=  1 


(41) 


where 


n 

o(t) = n  sin(27rrjt)  (42) 

i=  1 

and  the  interpolating  functions  FQ  (t)  at  the  origin  are  a  linear  combination 
of  G(f)/fJ,  j  =  1, . . .  ,  n,  chosen  so  that  Kf~1^  (0)  =  <5^/,  i,l  =  1, . . .  ,  n. 

4.3  Interpolators  at  Origin 

What  remains  for  this  interpolating  function  in  equation  (41)  is  to  find  ap¬ 
propriate  interpolators  at  the  origin.  To  remove  this  ambiguity,  I  construct 
interpolators  here.  Under  the  assumption  that  there  are  two  interpolators, 
then  G(z)  =  sin(27rz)  •  sin(27rnz),  where  a  is  the  square  root  of  any  prime 
number.  What  is  desired,  as  stated  in  the  theorem  in  section  4.2.2,  is  for 
the  interpolators  to  satisfy  the  following  condition:  {Ki(t)}  is  a  linear  com¬ 
bination  of  G{t)/tl,  i  =  1 ,n,  chosen  so  that  0)  =  <5y,  where 

j  =  !»•••  ,n- 

4.3.1  Two  Interpolators 

Assume  that  there  are  two  interpolators.  As  stated  earlier  in  equation  (42), 

G[t)  =  sin(27rt)  •  sin(2-7rat)  . 
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To  find  the  appropriate  interpolators  at  the  origin,  I  need  to  manipulate  this 
function.  Remember  that  the  Taylor  expansion  formula  of  a  sine  function  is 


sin(f)  =  t  -  +  0(t5) 

Thus,  inserting  equation  (43)  into  G(t)  obtains 

G(t)  =  sin(27rt)  •  sin(27ro;f) 

(27rt)3 


(43) 


27T  t 


3! 


+  0(t5) 


.  (27mt)3  n.5, 

27rat  — — — ■ h  0(t  ) 


=  47T2Crf2  — 


167T4a'3  167r4a 

+ 


3! 


3! 


3! 

t4  +  0(f6)  . 


And  so,  the  following  are  representations  for  the  linear  combinations  of 

G(t)/tl,  i  =  1. . . .  ,  n: 


Gify 

t 

Gjt) 

t 2 


47 x2at 


-  47 ra  — 


167T4Q:3  167T4Q: 

+ 


3! 

167T4Q:3 

3! 


+ 


3! 

167r4a 


3! 


t3  +  0(t5) 
t2  +  0(t4)  . 


Their  respective  derivatives  follow: 
d  G(t) 


dt  t 

d  G(t) 
dt  t2 


=  4:TTza  —  3 


167r4n3  167r4a 

+ 


3! 


=  2 


167r4a3  167r4a 

+ 


3! 


3! 


3! 

j 

t  +  0(t3) 


t 2  +  <D(t4) 


Next,  note  the  limits  of  the  functions  and  their  derivatives  at  t  =  0  — 


lim 

Git ) 

=  0 

t 

lim 

G(t) 

=  47r2a 

t- >0 

t2 

lim  4 

G(t) 

II 

t~*o  dt 

t 

lim  4 

G(t ) 

=  0  . 

t->o  dt 

t 2 

The  first  interpolant  K\  is  desired  to  be  such  that  Ki(0)  =  1  and  K[(0)  =  0. 
Also,  K2  is  desired  such  that  A^(0)  —  0  and  K'2{ 0)  =  1.  Setting  the  follow¬ 
ing  relations  satisfy  these  requirements: 


Ki  (t) 


GA0)f2 
2!  1 


G"  (0)  , 


(44) 


where,  conveniently,  one  can  see  that  G"(0)/2!  =  47r2a. 
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4.3.2  Arbitrary  Number  of  Interpolators 


It  is  convenient  to  further  construct  this  idea  to  an  arbitrary  multiple  of  in¬ 
terpolators.  In  this  scenario,  assume  that  there  are  n  interpolators  involved. 
Thus,  as  the  theorem  in  section  4.2.2  dictates,  there  are  n  sampling  rates 
Ui}  =  {v/po  =  1,  v'Pi,  •  •  •  ,  V^-T}'  where  ph  . . .  ,pn- 1  are  the  first  n  -  1 
primes,  which  generate  the  following  reconstruction  interpolant: 

n 

G(t)  =  JJsin(27iTjf)  .  (45) 

i=l 

Now,  the  process  of  finding  the  interpolants  at  the  origin  continues  as  be¬ 
fore.  With  the  use  of  the  Taylor  expansion  formula  for  the  sine  function, 
then 


G(t)  = 


n 

z-1  L 

n 

n 

i= l  L 


87r3r? 

27T  Tit - +  0(tb) 


3! 


2  nr. 


it_8Ty1.r.t3+0((5) 


3! 


=  2”^”  ( n  tn  -  2n+23f+2  ( n ^  tn+ 2 + °(*n+4) 


=  2  7T  y/piP^'-Pnt 


3! 

2^+2 ^n+2 


(46) 


3! 


\jvnn  •  •  -Pn 


f  n— 1  \ 

!  +  I Zpi)tn+2  +  °(tn+4) 

v  i= 1  / 


Now,  the  interpolants  at  the  origin  {AT*}  again  will  be  linear  combinations 
of  G(t)/t\  i  =  1, . . .  ,  n,  chosen  appropriately  so  that  (0)  =  Shj  ,where 

j  =  1, . . .  ,  n.  The  general  form  of  G(t)/tl  is 


Gify 

tk 


71  —  1 


nn+2  n+2  ( 

=  2VVWK-— -  J  ,/V:I»  --Vn  I  1  +  +  0(i”-t+i)  .  (47) 


«?(<) 


Let  jfffc(f)  =  G(t)/tk.  Then,  the  formula  for  the  jth  derivative  of  is 


2V,  (n  .  k)\,^pm---Pntn-k-3 


(n  —  k  —  j)\ 


2nT+\n(-I-)+2)!^TT^(1  +  g«)  t”_‘_2+2  +  Of*"'1-'"4)  •  <48> 


Allowing  t  to  approach  0,  the  following  results  are  obtained: 

2n7TnV£ 

G^(0) 


\im  Hk(t)  =  2n7Tn ^/'PlP2  ■  ■  ■  PnSnM 


n\ 


(49) 
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and 


lim  H{kj\t) 

t-»  o  fc 


Now,  using  these  evaluations  in  equation  (50)  as  before,  I  can  construct  a 
formula  of  appropriate  interpolants  using  the  set  {Hk}-  First,  note  that  each 
interpolant  is  of  the  form 


inn  (n  ~  fc)! 


2n7T 


{n-k-  j)\ 
2n+27rn+2  (n 


s/plP'l  '  '  ’  Priori, k-\-j 

—  k  +  2)! 


(50) 


3!  (n  —  k  —  j  +  2)! 


s/PlP~2  ■ ■ ■  Pn 


—  2 


k 

Kk(t )  =  Y,CiHi(t) 

2=1 

for  all  k  =  1, . . .  ,  n.  Note  that  the  c/s  vary  with  the  fc's.  As  required  of  the 
sampling  theorem,  coefficients  c,  are  needed  such  that  K^(0)  =  5k-i,j-  by 
linearity  of  differentiation,  I  have 


Kj*\0)  = 


2=1 

k 


C?; 


2=1 


f  77  —  7  V 

2n7r”  ~  T . ilTTVPlP2---pA,i+j 


2n+27rn+2  (n  —  i  +  2)!  ; ^ 

- 3!  +  5n'i+>~2 


=  E 


2=1 


r?  (n  — i)!  , - 


n— 1 


2n+27rra+2  (n  _  i)!  - /  i 

-Ci+2 - ^ ^  1  +  2_^Pi 

\  i= 1  / 


3!  ( n-i-j)! 


2=1 


(n-i-j)! 


$n,i+j 


n— 1 


=  2nnn^plp2  ■  ■  -Pv  •  j! 


227T2  j 

ci  -  Cj+2  oj  I  1  +  /  ,  Pi 
■  V  i=l  / 
2_2  /  "-1 


■>n,i+j 


Cn—j  Cfi—j+2' 


2 


3! 


!  +  5ZPi 


2=1 


G(n)(0)  . 


(n)! 

:  $k—\,j  • 


•j! 


^n— j  cn-j+2" 


22tt2 


3! 


2=1 


(51) 


In  equation  (51),  c^+i,  c*.+2  =  0  is  defined  so  the  summation  evaluates  cor¬ 
rectly.  Now,  remember  the  following  relation  must  hold:  K j,k  1)  (0)  =  1, 
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since 


Kri>(  o)  = 


(n)! 


(fc-1)! 


Cn— fc+1  cn—k+ 3" 


227T2 


3! 


'  n— 1  N 

1 + Ylpi 

\  i= 1  / 


•  (52) 


Since  (0)  =  0  if  j  k  —  1,  then 


cn—j  Cn—j+ 2" 


2V 

3! 


^  71  —  1  N 

V  2=1  / 


-0  . 


(53) 


Thus,  via  cancelation,  Cn-k+i  =  (G^ (0) / (n)\  ■  (k  —  l)!)-1.  Note,  also  from 
the  summations  just  mentioned,  that  there  is  the  restriction  n  —  j  <  k,  or 
restated  as  j  >  n  —  k.  So,  continuing  in  this  direction,  I  find  the  following 
recursive  equations: 


C/i— AH-l+2ra 


i^irq+EfaiVr 


G(")(0) 

(n)! 


‘  (k  —  1)! 


Cn— k+2+2m  —  0  . 


(54) 

(55) 


This  generates  interpolants  at  the  origin  with  the  desired  properties.  Thus, 
the  interpolants  can  be  formulated  as 


n—k 


2=1 


Ei£ot)/2)  (y(i  +  EfaiVi jj’V-Gffl 


G(")(  0) 
(n)! 


(, k  -  1)!  •  tn~k+1 


(56) 


Now,  with  these  interpolants,  I  can  use  this  new  reconstruction  method 
to  recover  a  function  of  a  larger  bandwidth  than  our  individual  sampling 
rates. 


4.4  Simulations 


Here  I  present  a  number  of  simulations  of  the  reconstructions  of  a  signal 
using  the  multirate  sampling  formula.  The  same  function  (signal)  used  in 
section  3.3,  f{t)  =  sin(27r£)/7rt.  Recall  that  this  is  an  Q  =  1  band-limited 
function.  The  classical  sampling  theory  stipulates  that  the  sampling  rate  be 
T  <  1/2 Q.  =  1/2.  The  classic  theory  dictates  that  the  Nyquist  rate  condi¬ 
tion  holds  to  successfully  reconstruct  the  function.  As  shown  earlier,  it  was 
sufficient  for  the  signal  to  be  sampled  at  Nyquist.  Now,  I  wish  to  sample 
slower  than  the  Nyquist  rate  and  still  reconstruct  the  function.  First,  the 
signal  will  be  sampled  twice  below  Nyquist  at  the  sampling  values 


/  1 

,  r2  = 

f  ±k  i 

|2  1  r-  1 

1  o 

(HV2)  J 

l/(lW5)  J 
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Note  that  these  sample  rates  are  both  slower  than  the  Nyquist  rate.  Also, 
note  that  these  are  appropriately  chosen,  i.e.,  1/(1  +  y/2)  +  \/2/(l  +  \/2)  =  1, 
and  the  ratio  of  the  two  samples  is  \/2  €  P.  Recall,  that  P  is  the  class  of 
irrationals  poorly  approximated  by  rationals.  Thus,  by  either  theorem  in 
section  4.2.2  or  4.2.3,  the  signal  can  be  reconstructed  by 

s(<)  =  t  /(x)G,(J((f_A,)  +  L/mtfiM + mm)i , 


with 

G(t) 

K\{t) 


K2(t) 


Essentially,  G(t)  looks  exactly  as  it  should,  i.e,.  the  product  of  two  sine  func¬ 
tions.  The  reconstruction  of  the  signal  is  in  figure  (4).  It  gives  the  appear¬ 
ance  of  completely  reconstructing  the  function  signal;  however,  upon  closer 
analysis,  minor  variations  can  be  seen  between  the  orignal  signal  and  the 
reconstruction.  These  points  of  aliasing  are  due  to  locations  in  the  sampling 
grid  where  the  samples  are  close.  Methods  of  eliminating  this  error  are  dis¬ 
cussed  by  Casey  and  Walnut  [12,19].  The  best  way  to  improve  this  sampling 
scheme  is  to  back  off  the  bandwidth  slightly,  i.e.,  increase  the  sampling  rates 
by  about  10  percent  or  more.  Then,  the  signal  can  be  reconstructed,  practi¬ 
cally,  within  an  acceptable  margin  of  error. 


=  sin 


2nt 


(1  +  \/2) 


sin 


2nV2t 

(1  +  V2) 


G{t) 

QG 01,2 
2!  1 

G(t) 

G"(0 ),  ' 


Figure  4.  Original  signal 
f(t)  and  bisample 
reconstruction  g(t). 


Time  (x) 
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Again,  as  in  the  classical  formula,  the  rate  is  critically  important  in  achiev¬ 
ing  an  accurate  reconstruction  of  the  signal.  Research  on  the  bounds  of 
these  rates  is  ongoing. 

Two  items  of  interest  are  the  interpolator  and  its  derivative  and,  most  im¬ 
portantly,  the  reconstruction  of  the  signal.  The  interpolator  and  its  deriva¬ 
tive  are  displayed  in  figure  5. 

In  figure  6, 1  demonstrate  the  reconstruction  of  the  same  signal  but  use  three 
noncommensurate  sampling  rates  below  Nyquist.  Here,  I  use  the  following 
sampling  rates: 


Ti 


±k 


(l+\/2+\/3) 


.  r2  = 


ifc 


V2 


,  r3  = 


±k 


(l+\/2+\/3) 


'(l+v^WS) 


Figure  6.  Original  signal 
f(t)  and  trisample 
reconstruction  g(t). 


Time  ( x ) 
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Note  that  these  rates  are  still  below  Nyquist.  The  reconstruction  formula 
used  is  the  same  as  presented  in  section  4.2.2.  The  only  real  point  of  interest 
to  present  is  the  interpolators  at  the  origin.  From  the  derivations  in  section 
4.3.2,  the  set  of  interpolators  is  shown  to  be 


m)  = 


m)  = 

K3(t)  = 


G(t)  +  •  (1  +  y/2  +  \/3)  ■  t2  •  G(t) 

G(3>(0),3 
3!  L 

G(t) 

G(3)(0)f9 

3!  1 

G(t) 

G(3)(0h  ' 

3!  1 


A  number  of  other  simulations  were  run  with  an  encouraging  success  rate, 
demonstrating  the  viability  of  this  multisampling  technique.  The  success 
rate  is  somewhat  improved  when  the  sampling  rates  are  more  equally 
spaced.  A  brief  discussion  of  applications  follows  in  the  next  section. 


4.5  Applications 


A  benefit  to  using  this  technique  is  that  the  individual  sampling  rates  are 
lower  than  the  traditional  Nyquist  sampling  rate  in  the  classical  sampling 
theorem.  For  example,  for  a  (1  +  aj-band-limited  function,  the  Nyquist  rate 
is  1  /2(l+n),  but  when  sampling  with  two  lattices,  the  sampling  rates  are 
only  1/2  and  l/2a.  Both  these  results  are  below  Nyquist.  The  generaliza¬ 
tion  shows  that  the  comparisons  between  the  sampling  rates  and  the  fre¬ 
quency  band  limit  are  more  outstanding  as  the  number  of  lattices  increase. 

Suppose  that  I  have  an  ft  band-limited  signal,  but  the  sampling  rate  is  re¬ 
stricted  to  1/2/?,  then  choose  the  first  n  primes  so  that 

(1  +Pi  +P2  H - Vpn)'—  >ft-  (57) 

Pn 

Then,  n  +  1  sampling  lattices  can  be  constructed  via  the  theorem  in  section 
4.2.3  that  can  effectively  reconstruct  /  with  sampling  rates  no  faster  than 
1/2/?. 

Unfortunately,  this  structure  has  some  difficulties.  The  first  is  that  the  sam¬ 
pling  grid  is  extremely  rigid;  perturbations  on  the  lattices  result  in  losing 
information.  Also,  because  the  sampling  points  might  become  arbitrarily 
close,  ripples  occur  at  points  near  each  other.  If  the  bandwidth  is  backed 
off,  the  simulations  results  improve,  but  exact  bounds  are  unknown. 

These  results  also  extend  to  higher  dimensions,  where  one  can  think  of 
the  grids  and  the  interpolating  function  as  Cartesian  products.  Casey  and 
Sadler  [20]  have  used  these  techniques  to  develop  new  analog-digital  trans¬ 
formers  for  signal  processing.  In  their  signal  processing  literature,  these  are 
known  as  A-D  converters. 
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Appendix  A.  Proof  of  the  Polynomial  Expansion 


I  shall  show  the  result  in  section  3.4.2  of  the  main  report  via  induction  on 
the  degree  of  the  polynomial. 

Proof.  Suppose  p(Q  is  a  one  degree  polynomial  with  leading  coefficient  1 
and  one  simple  zero,  then  p(()  =  (  —  r\.  Thus,  clearly, 

~7n  =  ■  (A-D 

p(Q  C  -  n 

Now,  assuming  that  this  statement  is  true  for  any  k  degree  polynomial  with 
leading  coefficient  1  and  k  distinct  simple  zeros,  where  k  =  1 , . . .  ,  s,  I  want 
to  show  the  same  result  for  an  s  +  1  degree  polynomial  with  s  simple  zeros. 
Suppose  p(C)  is  such  a  polynomial.  Then,  I  have  a  representation  for  p  using 
the  set  of  its  zeros: 

i>(0  =  ntc-rj)  .  (A-2) 

3=  0 

I  can  equivalently  state  this  as  p( C)  =  (C  —  J’o)  •  (l{z),  where  q  is  an  ,s  de¬ 
gree  polynomial  with  only  s  simple  zeros.  Thus,  I  have  a  formula  by  the 
assumption  in  the  induction  for  q.  Using  this,  I  obtain  the  following: 

1  _  1  1 

p(0  (C-ro)  q(z) 


(C  -  ro)  (<  -  Tj)  ■  q'irf) 

s  1 

_  ^  ( C  -  ro)(C  -  U'MU')  ' 


31 


Now,  let  m(C)  =  «  -  r0)(C  ~  r^q'irj),  then  m'( ()  =  ((  -  rj)q'(rj )  +  (C~ 
ro)q'(rj)-  Also,  note  that  p'{Q  =  q(()  +  (C  -  r0)q'(Q.  So  continuing,  I  have 


1 

p(  0 


=  y _ _ _ 

y  (C-  r0)(C  -  rj)q'(rj) 


-£ 


(C  -  r0)m'(r0)  (C  - 


_  A _ 1 _ + _ 1 _ 

~  “  (C  -  ^o)(ro  -  rj)q'(rj)  (C  -  rj)(rj  -  r0)q'{r3) 

1  V— ^  1  g  ^  J 

(C  -  r0)  ^  (ro  -  rjWirj)  +  ^  (C  - 

-  1  1  A  1 

(C  -  ro)  q(ro)  +  y  (C  -  rj)v'{rj) 


(A-4) 


1 

(C  -  r0) 


i  A  l 
p'(r 0)  (c  -  ' 


Equation  A-4  gives  the  desired  result. 


Appendix  B.  MATLAB  Code 


The  following  MATLAB  *.m  files  were  used  to  perform  the  simulations.  I 
used  MATLAB  version  5.3  by  The  Mathworks,  Inc. 


B-l  Shannon.m 

'/SHANNON  EXAMPLE  '/,  This  m-file  generates  plots  for  the  demonstration  of 
Shannon’s  formula.  Three  plots  will  be  created.  The  first  will  show  the  signal  and  its 
Fourier  transform,  the  second  will  demonstrate  the  reconstruction  procedure,  and 
the  last  plot  will  show  the  effects  of  aliasing  in  the  sampling. 

'/,  Terrence  Hoore,  3/6/2000 

'/,  Revision:  1.3  Date:  7/5/2000  18:59:33 

'/.begin  of  code 
'/generate  the  domain 
x  =  linspace(-5, 5,1000); 

'/.the  signal  function 
y  =  sin(2*pi*x)./(pi*x); 

'/.plot  of  signal 
figure; 

subplot(3,2,3); 

plot(x,y); 

title(  'original  signal  f(x)’); 

xlabel(’x’); 

ylabel(,f(x)0; 

'/.the  transform  of  the  signal 
xl  =  [-2  -1-1112]; 
yl  =  [0  0  1  1  0  0] ; 

'/.plot  of  the  transform 
subplot(3,2,4); 
plot(xl,yl); 
axis([-2  2-12]); 

title( Courier  transform  f(\  omega)’); 
xlabel(’\  omega’); 
ylabel(’f(\  omega)’); 
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lour  re contraction  formula  only  for  this  signal 

'/.omega  =  1 

s=l; 

z=0; 

count=0; 
for  n=-20:20 
if  n=0 

z^z+sintpi+n/s) . /(pi*n/(2*s))  .*sin(pi*(2*si«:x-n))  ./(pi*(2*s*x-n)); 
count=count+l; 
end 
end 

z=z+2*sin(pi*(2*s*x)) ./(pi*(2^s*x)) ; 
figure; 

plot(x)y,,-,,xJz,,:}); 

hold; 

'/.the  sampled  data 
N  =  -20:20; 

N  =  N/(2*s); 

yN  =  sin(2*pi*N)./(pi*lt); 
plot(N,yN,  V);plot(0}2,  V); 
axis([-5  5  -0.5  2.5]) 
xlabel(’x’); 

legend( ’f (x) ’ , ’g(x) ’ , ’ samples’ ) ; 

'/.omega  =  0.8  (under sampling  -  aliasing) 

s=0.8; 

z=0; 

count=0; 
for  n=-20:20 
if  n  =0 

z=z+sin(pi*n/s) . /(pi*n/(2*s)) . *sin(pi*(2*s*x-n)) ./(pi*(2*s*x-n)) ; 
count=count+l; 
end 
end 

z=z+2*sin(pi*(2*s*x) ) .  /(pi*(2*s*x)) ; 
figure; 

plot(x,y,’-’,x,z,’:’); 

hold; 

'/.the  sampled  data 
N  =  -20:20; 

N  =  N/(2*s); 

yN  =  sin(2*pi*N)./(pi*N); 
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plot(NjVo,);plot(0,2,V); 
axis ([-5  5  -1  2.5]) 
xlabel(Y); 

legend(  ’f  (x) ’ ,  Jg(x) * ,  Samples  ’) ; 

'/.omega  =  1.2  (oversampling  --  aliasing) 

8=1.2; 

z=0; 

count=0 ; 
for  n=-20:20 
if  n  =0 

z=z+sin(pi*n/s)./(pi*n/(2*s)).*sin(pi^2*s*x-n))./(pi*(2*s*x-n)); 

count=count+l; 

end 

end 

z=z<!*sin(pi*(2*sn))./(pi*(2*s*x)); 

figure; 

plot(x,y,,-,,xJz/:  ’); 
hold; 

'/.the  sampled  date 
H  =  -20:20; 

N  =  N/(2*s); 

yN  =  sin(2*pi+N) . / (pi*N) ; 
plot(N,yN,  V);plot(0,2,  V); 
axis ([-5  5  *0.5  2.5]) 

title(,reconstruction  function  g(x)  and  signal  f(xV); 
xlabelf'x1) ; 
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B-2  MultirateShannon.m 

t  =  linspace(-5,5, 1000) ; 

y  =  sin(2*pi*t) ./(pi*t); 

rO  =  1  +  sqrt(2); 
r  =  [1  sqrt(2)]; 
r  =  r./rO; 

N  =  20; 

rates  =  l./(2*r); 

Ns  =  round (K* ( 1 . /rates)  +1); 
tsamp  =  zeros(length(r)  ,2jtcmax(Ks)+l); 
for  i  =  1 : length(r) 

tsamp(i,:)  =  [(-Ms(i):Ns(i))  zeros(l,2*max(Ns)-2*Ns(i))] ; 
end 

samp  =  diag(rates)*tsamp; 
samp  =  samp( : ) ; 
indz  =  find(samp  ==  0); 
samp(indz)  =  []; 
lambda  =  sort([samp. *]); 
lambda  =  lambda;*/,  +  .001; 

ST  =  sin(2*pi*r. J*t); 

St  =  prod(ST,l); 

CT  =  2*pi^diag(r)*cos(2*pi*r. "*t) ; 

CT  =  flipud(CT); 

Sp  =  sum(ST. +CT, 1) ; 

figure; 

subplot(2,l,l); 

plot (t, St); 
subplot(2,l,2); 
plot (t ,Sp) ; 

SI  =  sin(2*pi*r. ’♦lambda); 

Cl  =  2*pi*diag(r)*cos(2*pi*r . '^lambda) ; 

Cl  =  flipud(Cl); 
bp  =  SI . *C1 ; 
bp  =  sum(bp,l); 

'/.sample  values  to  be  interpolated 
tp  =  sin(2*pi*lambda)7(pi*lainbda); 
fp  =  tp./bp; 
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temp  =  zeros(size(t)); 

for  i  =  1: length (lambda) 

temp  =  temp  +  fp(i) . /(t-lambda(i)) ; 

end 

f  =  temp^St; 

ST2  =  2*4*pi*pi*prod(r,2); 

K1  =  St./(t.*t); 

K1  =  K1/ST2*2; 

K2  =  St./(t); 

K2  =  K2/ST2+2; 


f  =  f  -  2*K1  +  0*K2; 
figure; 

plot(t>y,,-,,t,f  ,’:’); 
hold; 

'/.the  sampled  data 
N  =  lambda; 

yN  =  sin(2*pi*N). /(pi*N); 
plot(N,yM,  V);plot(0,2,  V); 
axis([-5  5  >1  2.5]) 
xlabel(’x’); 

legend( ?f (x) ’ , ’g(x) * , ’samples’) ; 
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